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q'; abstract 

Shock waves developed during the formation and evolution of cosmic structures 
, are key features encoding crucial information on the hierarchical formation of the 

04 Universe. We present the analysis of an Eulerian adaptive mesh refinement (AMR) 

' I hydrodynamical and N-body simulation in a ACDM cosmology especially focused on 

^ . the study of cosmological shock waves. The combination of a shock-capturing algo- 

rithm together with the use of a halo finder allows us to study the morphological 
\ structures of the shock patterns, the statistical properties of shocked cells, and the 

correlations between the cosmological shock waves appearing at different scales and 
the properties of the haloes harbouring them. According to their localisation with 
respect to the population of haloes in the simulation, shocks can be split into two 
broad classes: internal weak shocks related with evolutionary events within haloes, 
I and external strong shocks associated with large-scale events. The shocks segrega- 

. tion according with their characteristic sizes is also visible in the shock distribution 

function. This function contains information on the abundances and strength of the 
different shocks, and it can be fitted by a double power law with a break in the slope 
around a Mach number of 20. We introduce a generalised scaling relation that cor- 
relates the average Mach numbers within the virial radius of haloes and their virial 
masses. In this plane, Mach number - virial mass, two well-differentiated regimes ap- 
pear. Haloes occupy different areas of such plane according to their early evolutionary 
histories: those haloes with a relatively quiet evolution have an almost constant Mach 
number independently of their masses, whereas haloes undergoing significant merger 
events very early in their evolution show a linear dependence with their masses. At 
\ high redshift, the distribution of haloes in this plane forms an L-like pattern that 

evolves with time bending the vertical branch towards the horizontal one. We prove 
that this behaviour is produced by haloes reaching low average Mach numbers as they 
evolve towards a virialised state. The analysis of the propagation speed and size of 
the shock waves developed around haloes could give some hints on the formation time 
and main features of the haloes. The possible future detection of all sort of cosmic 
shocks by the forthcoming telescopes, such as the SKA, could open a new window to 
indirectly unravel some of the main features of haloes and, therefore, the formation 
process of the cosmic structures. 

Key words: hydrodynamics - methods: numerical - galaxy clusters - large-scale 
structure of Universe - shock waves - Cosmology 
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1 INTRODUCTION 



The cosmological shock waves develop as a consequence of 
the hierarchical formation of structures in the Universe and, 
* e-mail: susana.planellesOoats.inaf.it therefore, they are crucial ingredients in a unified picture of 

t c-mail: vicent. quilisOuv.es the formation of cosmological structures. In a first phase, 
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gravitational energy associated to the collapse of dark mat- 
ter haloes corresponding to galaxy clusters and galaxies is 
transformed into internal energy of the intra cluster medium 
(ICM) and inter galactic medium (IGM) gaseous compo- 
nents. In the following phases, the evolution of those struc- 
tures also produces mergers and accretion phenomena that 
modify the energetic balance of the gas through shock waves. 
Therefore, the shocks associated to cosmic structure forma- 
tion and evolution encode information about the formation 
of the structures and their thermal impact on the gas. 

The cosmological shocks can be broadly classified into 
two categories (e.g., IRvu et al] l2003h : external and inter- 
nal. External shocks surround filaments, sheets, and haloes, 
while internal shocks are located within the regions bound 
by external shocks and are created by flow motions cor- 
related with structure formation and evolution. On large 
scales, the thermal history of galaxy clusters is dominated 
by the infall of material onto dark matter haloes and the 
conversion of gravitational energy into thermal energy of the 
gas. This process occurs through the heating of the gas via 
strong (external) accret ion shock s surrounding galaxy clus- 
ters and filaments (e.g .,| Rvu et al. 2003: Mi niati ct al.,, 200ll : 
IPfrommer et al.ll2006l ). Inside collapsed structures, weaker 
(internal) shocks can be subdivided in three different classes: 
(i) accretion shocks caused by infalling gas in cosmic struc- 
tures, (ii) merger shocks resulting from merging haloes, and 
(iii) random flow shocks inside nonlinear structures pro- 
duced during hierarchical clustering. These internal shocks 
contribute to the virialization of haloes. 

The role of shocks in cosmological structures has been 
studied from different complementary approaches. From 
the observational point of view, strong shocks usually de- 
velop in the external low-density regions of galaxy clus- 
ters where observable emission, such as X-ray emission, is 
weak. As a consequence, from an observational point of 
view, detecting shocks on large scales is still very chal- 
lenging. Despite these difficulties, few large cosmic shocks 
have been positiv ely detected by me ans of the so-called 
radio relics (e.g., lEnsslin et al.l Il998l ). commonly defined 
as elongated radio emission not associated with the clus- 
ter centre or an active cluster radio galaxy. These radio 
relics are found in ~ 30 clusters (e.g., Bagchi et al.l 20061: 
Venturi et al.ll2007l : ISonafede et aLll2009r i van Weeren et al.l 



20091 ). Peripheral radio relics on the outskirts of some 



massive clusters are the brightest of the cluster shock 
structures and are interpreted as externally propagating 
merger shocks that reach higher Mach number^ as they 
enter the low-density regions of the ICM. Concerning in- 
ternal shocks, in a few cases, shocks driven by merging 
events have been observed with very low (~ 1.5 — 3) 
Mach numbers (e.g., Markevitch. Sarazin fc Vikhlininlll999l: 



iMarkevitch et al.ll2002l : iMarkevitch fc Vikhlininll200it ). 

From the theoretical point of view, several at- 
tempts to st udy shocks were d o ne us in g semi-analytical 
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Gabici fc Bias 



Pavlidou fc Fields! 
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appr oa ches (|Press fc Schechtej 1 1974 \ Sheth fc TormenI 



120061 : iFuiita fc SarazinI l200lF 
However, the complexity of the 



considered scenario proved the use of numerical approaches 



^ The Mach number, which characterises the strength of shocks, 
is defined in Eq. [T] 



as the most suited to study and characterise the shocks 
produced during the formation and evolution of cosmic 
structures. There have been numeric al studies of shocks us- 
ing b o th Eulerian, "s i ngle-g r id" (e.g., Quilis. Ibawez fc Saez 



19981 : iMiniati et all l200ol: IRvu et al 



20071 : IVazza. Brunetti fc Ghelleij |2009|) and "AM R- 



20031 : iKang et al 



^^o ^ches Jskillma; r:n!n ^t ^ (vazza et al. l 



grid; 
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2010h. as well as SPH codes (e.g., JPfrommer et ;a.l.l 120061 : 
Pfrommer. Enfilin fc Springell I2OO8I : iHoeft et all booj ). 
Even when the debate about the benefits and draw- 
backs of different numerical techniques is still open (e.g., 
lAgertz et a"l]|2007f ). the inherent shock-capturing properties 
of the Eulerian techniques based on Riemann solvers, able 
to deal with shocks by construction, seem to position this 
kind of tools as the recommended o nes when t a ckling 
shock related processes. In this regard, IVazza et all (|201ll ) 
presented a numerical study of non-radiative cosmological 
simulations, at various resolutions and performed with 
different cosmological codes, comparing the properties of 
thermal gas and shock waves in large scale structures. 
Even when the bulk of thermal and shock properties were 
reasonably in agreement between the Eulerian and the 
SPH codes, they also reported some significant differences 
between them like, for instance, differences of large factors 
(rsj 10 — 100) in the values of average Mach numbers and 
shock thermal energy flux in the most rarefied regions of 
the simulations, significantly different phase diagrams of 
shocked cells in grid codes compared to SPH, or sizable 
differences in the morphologies of accretion shocks between 
grid and SPH methods. 

Pioneering attempts to characterise shock waves 

in cosmological simulation s we re carried out by 

iQuilis. Ibanez fc Saej (|l998l ) and iMiniati et all (|2000l ). 
They employed fix grid Eulerian simulations and shock- 
detecting schemes based on jumps in the main thermo- 
dynamical quantities. Later works adopted more refined 
shock-detecting algorithms and were more focused onto 
the dis t ributi o n of energy diss i pated at shocks (e.g. , 
Miniatil [20o3: IRvu et al l l2003l : IPfrommer et all I2OO6I : 



Vazza. Brunetti fc GhelieiT |2009| ). Despite the important 



advances, in the first works using uniform grid-based codes, 
it was not possible to cover the spatial resolutions required 
to describe both the complex flows within haloes and their 
coupling to large-sc ale structures. Further improvement 
has been reached by ISkillman et al.l (|2008l ) using a shock- 
detecting scheme looking for shocks in the direction of the 
temperature gradients on an AMR grid. 

These numerical simulations have begun to re- 
veal a rich net work of sho c k stru c tures thro u ghout 
the ICM (e.g.. IMiniati et all [2001bl: IRvu et all 20031 : 



Pfrommer. Enfihn fc Springell I2OO8I: iBattaglia et all bood : 



Skillman et al.ll2010l : IVazza et al.ll2010l ). However, in spite 
of all previous works, the identification and characterisation 
of shocks is still challenging due to the complex dynamics in- 
volved in the formation and evolution of cosmological struc- 
tures and to the large dynamical range needed to describe 
all the scales involved by shocks. 

Our purpose in the present paper is to pursue the anal- 
yses of the main properties of the shock waves developed 
during the evolution of a high resolution hydrodynamical 
and N-body simulation of a large cosmological volume per- 
formed with an AMR cosmological code. In addition, we will 
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put especial emphasis in analysing the existing connection 
between the cosmological shock waves and the population 
of haloes. In order to do so, we have developed a numerical 
algorithm able of detecting and characterising shocks in 3-D 
AMR simulations. The use of AMR hydro codes turns out 
to be crucial so as to obtain good dynamical ranges with an 
advanced hydrodynamical algorithm that is able to capture 
shocks very accurately. 

The paper is organised as follows. In Section 2, we 
present the technical details describing the simulation and 
the shock-finding algorithm. The results of the analysis of 
the simulation are presented in Section 3. Finally, in Section 
4, we summarise and discuss our results. 



2 DETECTING SHOCK WAVES 
2.1 Simulation details 

The simulation used in this p aper was per formed with the 
cosmological code MASCLET (|Quilisll20oJ ). This code cou- 
ples an Eulerian approach based on the so-called hiqh - 
resolution shock- capturing techniques (e.g., lLeVeQu3ll992l ) 
for describing the gaseous component, with a multigrid par- 
ticle mesh N-body scheme for evolving the coUisionless com- 
ponents (dark matter and stars). Gas, dark matter, and stars 
are coupled by the gravity solver. Both schemes benefit of 
using an adaptive mesh refinement (AMR) strategy, which 
permits to gain spatial and temporal resolution. 

The numerical simulation was run assuming a spa- 
tially flat ACDM cosmology, with the following cosmo- 
logical parameters: matter density parameter, Jim = 0.25; 
cosmological constant, Ha = A/3H^ — 0.75; baryon den- 
sity parameter, ilb = 0.045; reduced Hubble constant, 
h = Ho/WOkms-^ Mpc-^ = 0.73; power spectrum index. 
Us = 1; and power spectrum normalisation, erg — 0.8. 

The initial conditions w ere set up at z = 50, u sing a 
CDM transfer function from lEisenstein fc Hul (|l998h . for a 
cube of comoving side length 64 Mpc discretised in 512^ cu- 
bical cells. 

A first level of refinement (level / = 1) for the AMR 
scheme was set up from the initial conditions by selecting 
regions satisfying certain refining criteria, when evolved - 
until non-linear phase- using the Zel'dovich approximation. 
The dark matter component in the initial refined region was 
sampled with dark matter particles eight times lighter than 
those used in regions covered only by the coarse grid (level 
1 = 0). During the evolution, regions on the different grids 
are refined based on the local baryonic and dark matter 
densities. Thus, a cell is refined, independently of the refine- 
ment level, if its dark matter (gaseous) mass is larger than 
4.74 X 10* M© (1.04 X 10* M©). This is equivalent to refine 
the cell if its density increases a factor of eight. The ratio 
between the cell sizes for a given level (/ + 1) and its parent 
level (l) is, in our AMR implementation, Axi+i/ Axi = 1/2. 
This is a compromise value between the gain in resolution 
and possible numerical instabilities. 

The simulation presented in this paper uses a maximum 
of six levels {I = 6) of refinement, which gives a peak physical 
spatial resolution of ~ 4 kpc. For the dark matter we consider 
two particles species which correspond to the particles on 
the coarse grid and the particles within the first level of 



refinement at the initial conditions. The best mass resolution 
is ~ 2.7 X 10* Mq, equivalent to distribute 512^ particles 
in the whole box. 

Our simulation includes cooling and heating processes 
which take into acco unt inverse Compton and free-free 
cooling, UV heating (|Haart fc Madaul Il996l l at z ~ 6, 
atomic and molecular cooling for a primordial gas, and 
star formation. In order to compute the abundances of 
each species [H, He, , He^, He^^), we assume that 
the gas is optically thin and i n ionization equilibrium, but 
not i n therma l equilibrium (|Katz. Weinberg fc Hernguistl 
1 19961 : [ Theuns et al. 199 j). The tabulated cooling rates were 
taken from Sutherland fc Dopital (| 19931 ) assuming a con- 
stant metallicity 0.3 relative to solar. The cooling curve was 
truncated below temperatures of 10* K. The cooling and 
heating were included in the energy equation (see Eq. 3 in 
Ifluilis 20()|) as extra source terms. 

The star formation is intr oduced in the M AS- 
CLET code following the ideas of lYepes et al.l (|l997D and 
ISpringel fc Hernguistl (|2003l ). In our particular implementa- 
tion, we assume that cold gas in a cell is transformed into 
star particles on a characteristic time scale according to 
p, = —f> = (1 — /?) p/t,{p) where p and p, are the gas and 
star densities, respectively. The parameter /3 stands for the 
mass fraction of massive stars (> 8Mq) that explode as su- 
pernovae, and therefore return to the gas component in the 
cells. We adopt /3 = 0.1, a value compatible with a Salpeter 
IMF. For the characteristic star formation time, we make 
the common as sumption ^t(p) = io{p/ Pth)^^^^ ' equivalent 
to p« = p^'^ /t* (|Kennicuttl 19981 ) . In this way, we introduce a 
dependence on the local dynamical time of the gas and two 
parameters, the density threshold for star formation (p^^) 
and the corresponding characteristic time scale {tg). In our 
simulation, we take tg = 2Gyr and p^^ = 2 x 10~^^ gcin~^. 
From the energetic point of view, we consider that each su- 
pernova dumps in the original cell 10^^ erg of thermal en- 
ergy. 

In the practical implementation, we assume that star 
formation occurs once every global time step, Ati^o, and 
only in the cells at the highest level of refinement. Those 
cells at this level of refinement, where the gas temperature 
drops below T < 2 x 10* K, the divergence of the velocity 
of the gas is V ■ v < 0, and the gas density is p > Pj,^ = 
2 X W~^^ g cm~^ , are suitable to form stars. In these cells, 
coUisionless star particles with mass m* = p^Ati^gAxf are 
formed. In order to avoid sudden changes in the gas density, 
an extra condition restricts the mass of the star particles to 
be m, = min(m,, ^rrigas), where mgas is the total gas mass 
in the considered cell. 

According with the previously described scheme for the 
star formation in our simulations, we include only one sim- 
ple case of feedback, that is, the contribution in mass and 
energy from type II supernovae. No feedback from active 
galactic nucleus (AGN) or type la supernovae are consid- 
ered. Although this point could be improved, how to model 
AGN feedback in c osmological simula tions is still a matter 
of discussion (e.g., iFabian et al.|[201oh . Besides, we do not 
expect a dramatic change in our results linked with the in- 
clusion or not of feedback, since the formation and evolution 
of the cosmological shock waves are supposed to be mainly 
driven by the assembly of cosmic structures and, therefore, 
gravity should be the main ingredient to model. 
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A final technical consideration concerning the numeri- 
cal resolution of the simulation must be done. We have not 
performed any resolution test and assume that our numeri- 
cal resolution is adequate to study the processes related with 
shocks. We base our assumption in the work bv lVazza et al.l 
()201lh . where authors compare the properties of shocks us- 
ing different cosmological codes. We must point out that our 
numerical resolution - in cell size and particles masses - is 
much better than in all the simulations considered in Vazza 
et al. (2011) (see Table 1 of this reference), where the au- 
thors conclude that a very good convergence, better than a 
^ 10 per cent level, in the most important shock statistics 
is expected for spatial resolutions of 50 — 100 kpc, which 
are above than our best resolution. 



2.2 Basic relations 

Shocks produce irreversible changes in the gas of the cos- 
mic structures. As a consequence, the formation of a shock 
in a cosmological volume produces a jump in all the ther- 
modynamical quantities. If we assume that the pre-shocked 
medium is at rest and in thermal and pressure equilibrium, 
the pre-shock and post-shock values for any of the hydro- 
dynamical variables are unambiguously related to the shock 
Mach number, M. The Mach number, which characterises 
the strength of a shock, is given by 



M = Vs/Cs 



(1) 



where Vs is the shock speed and Cs is the sound speed ahead 
of the shock. 

All the information needed to evaluate M is contained 
in the Rankine-Hugoniot jump conditions. If the adiabatic 
index is set to 7 = 5/3 we obtain, for the density (p), the 
temperature (T), and the entropy {S = T/p^ "^), the well- 
known relations fe.g.. iLandau fc Lifshitdll966f l: 
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(3) 



(4) 



16X2 

{5M^ ^1){M'' + 3) +3 2/3 

with indices 1, 2 referring to pre- and post-shock quantities, 
respectively. 

The Mach number can be obtained from the jumps in 
any of the hydrodynamical variables (Eqs. [2}[4]) or from a 
combination of them. In the case of relatively large Mach 
numbers, since the value of the density jump saturates at 
P2/P1 = 4 (Eq. [2|, strong shocks cannot be detected from 
density jumps. As a consequence, the most effective meth- 
ods to measure A4 are those considering temperature and 
entropy jumps. 

2.3 Shock- finding algoritlim 

Any shock-finding method relies on accurately identifying 
and quantifying the strength of shocks. Although, during 
the simulations using an Eulerian approach, the shocks are 
automatically detected by the Riemann solver within the hy- 
drodynamical routine, the analysis of these shocks requires 
additional considerations. 



Recent works on algorithms to find and track shocks 
in grid simulations have stablished two procedures. Both 
methods require a previous step to mark the cells which are 
involved in a shock at the time corresponding to the anal- 
ysed output. Once the shocked cells have been identified, 
the two meth ods differ in how to define the sho cks. The first 
method (e.g.. I Vazza. Brunetti fc Ghellejlioogl ') relies on fol- 
lowing temperature or velocity jumps across the previously 
tagged shocked cells along the different coordinate axes. We 
define this method as th e coordiante or directional-splittin g 
method. According with lVazza. Brunetti fc GheUeij (|2009h . 
the results seem to be quite insensitive to the use of tem- 
perature or velocity jumps, with minor differences in the 
most rarefied environments. This kind of methods has a mi- 
nor drawback when examining shocks that do not propagate 
along the coordinate axes. The second group of methods to 
find shocks use the local temperature gradie nt to figure out 
the d irection of the shock propagation (e.g.. ISkillman et al.l 
l2008h . These methods would be more precise when dealing 
with complex fiows or weak shocks. So as to minimise a pos- 
sible smearing of the shocks produced by the coordinate- 
splitting or other uncontrolled numerical effects, the so- 
called coordinate-splitting methods define the final Mach 
number as the average of the Mach numbers along each co- 
ordinate axis. With this implementation, the differences be- 
tween schemes based on a precise estimate of the direction of 
the shock propagation and others based on the directional- 
splitting are not expected to be significative. 

Our shock-finding algorithm is based on tempera- 
ture jumps and, therefore , comp arable to the method in 
IVazza. Brunetti fc GhelleJ (|2009l '). Thus, the main steps of 
the shock-finding algorithm used in this paper are the fol- 
lowing: 

(i) All cells within the computational volume are classi- 
fied as tentative shocked or not shocked. A cell is labelled 
as tentative shocked if the fiuid inside the cell meets the 
following requirements: 



V ■ V < 
VT ■ V5 > 



(5) 
(6) 



where v, T and S are, respectively, the velocity field, the 
temperature and the entropy of the fiuid (gas) within the 
cell. 

(ii) Among all the tentative shocked cells, the cell where 
V ■ V is minimum is flagged as the first shocked cell. 

(iii) The sense of shock propagation with respect to the 
computational grid is computed by moving outwards from 
the previously identified shocked cell along each of the three 
coordinate axes. This extension is done meanwhile there are 
tentative shocked cells and the temperature and the den- 
sity in these tentative shocked cells satisfy that T2 > Ti 
and p2 > pi- The subscripts 1 and 2 refer to the tenta- 
tive shocked cells ahead and behind of the central shock cell 
in the shock reference frame, respectively. The shock dis- 
continuities caught in the simulations by this method are 
typically spread over a few cells. 

(iv) Once the furthest shocked cell ahead and behind of 
the central cell of the shock - that we denote as the pre- 
and post-shock cells - in each coordinated direction are 
found, the temperatures Ti (pre-shock) and T2 (post-shock) 
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of these cells are taken and three different Mach numbers 
(one for each direction) are calculated from Eq. [S] 

(v) The strength of the shock is computed by com- 
bining the Mach numbers measured along the three co- 
ordinate axes: M = {Ml + Ml + MIY^'^, which min- 
imises projection effects in cas e of diagonal shocks (e.g., 
IVazza. Brunetti fc Ghelleill2009l '). 

(vi) The shocked cell drops out of the list of the tentative 
shocked cells, and the process is iteratively repeated focus- 
ing on the cell with minimum V ■ v among the remaining 
tentative shocked cells. 

The whole procedure identifies all the shocked cells 
within the computational box and obtains their Mach num- 
ber. The assembly of all these shocked cells defines the char- 
acteristic shock surfaces associated to shock waves. 

When this procedure is applied to an AMR simulation, 
the analysis is carried out in a hierarchical fashion, first on 
the most highly refined grids, moving down to progressively 
coarser levels of resolution. Given that this procedure is ap- 
plied, independently, at each level of refinement of the simu- 
lation, the algorithm is able to find, in a natural way, shock 
waves related to different spatial scales provided by the sim- 
ulation itself. Although this is a very useful and interesting 
feature of our approach, the use of an AMR grid makes 
the process more complicated compared with a fix grid ap- 
proach, especially when shocks extend over several subgrids 
at different levels. In this case, the algorithm goes on looking 
for the boundaries of the shock wave on the parent subgrid 
- in a lower level of refinement - that contains the original 
subgrid on which the shock was originally identified. This 
mechanism can be applied recursively on several levels if 
needed. 

In all the analysis performed in following sections of this 
paper, in order to avoid noisy shock patterns with very low 
Mach numbers, we have considered a Mach number mini- 
mum threshold equal to 1.3. 



3 RESULTS 

3.1 Distribution of large-scale shocks 

Shoc ks fill the simulated volume in a very complex pattern 
(e.g., iMiniati et al] I2OO0I : IRvu et all 120031 '). Figure [1] illus- 
trates typical structures found in large-scale cosmological 
simulations. The different two-dimensional panels (64 Mpc 
side length) represent the following integrated quantities 
along slices of 10 Mpc thickness at z = 0: the gas over- 
density (upper left), the gas temperature (upper right), the 
gas entropy (lower left) and the Mach number (lower right). 
All the panels show the logarithm of the different quanti- 
ties. The four projections are centred at the position of the 
largest cluster in the simulatior0 {Myir ~ 8.0 X 10^'' Me and 
Rjjir ~ 2.4 Mpc) which is almost at the centre of the box. 

The panel presenting the gas overdensity in Fig.[T]shows 
the common picture of galaxy clusters connected through 
filaments. The shock waves, if existing, are not appreciable 
in this picture. The panel displaying the gas temperature 

^ See Section 3.2 for further details on the halo population ob- 
tained in the present simulation. 



clearly shows high temperature regions, where the gas has 
been heated up by shock moving outwards from the centre 
of the structure, and low temperature spots associated to 
regions where the gas cooled down during the cosmic evolu- 
tion. The temperature plot gives a first preliminary idea of 
the regions which have been shocked. The entropy panel, as 
it would be expected, is a step forward highlighting the com- 
plex structure of shocks. In this particular panel, the bound- 
aries of high temperature regions appear as high entropy 
zones indicating the presence of a shock or strong gradient. 
For completeness, the last panel shows the Mach number as- 
sociated to the studied region. Now, the complex structure 
of shocks is clearly visible, showing a remarkable correlation 
with the entropy panel, although much sharper in this case. 
The Mach number map also allows to distinguish between 
high-Mach number shocks, wrapping the larger structures, 
and the low-Mach number shocks associated with internal 
small scale structures. 

The cosmic evolution of shock strengths encodes a valu- 
able information about the thermal history of the baryonic 
component of the Universe. In Fig. [21 the filling volume fac- 
tor of shocks - that is, the fraction of the computational vol- 
ume occupied by shocked cells - is plotted versus redshift 
evolution (left panel), whereas the right panel shows the 
volume-weighted average Mach number for the same time 
evolution. A warning notice must be done at this point be- 
cause, as already mentioned in Sec. 2.1, a resolution test 
has not been performed. Therefore, these results could be 
affected by resolution effects. 

Early evolution (z > 6) of the shocked volume evidences 
a non-negligible fraction of the computational box involved 
in shocks (~ 27%). This quantity increases towards its max- 
imum at 2 ~ 5 when a step decline of the filling factor of 
shocks starts. The explanation of such behaviour is asso- 
ciated with huge amounts of gas falling into the potential 
wells created by cosmic structures - at different scales - at 
the beginning of their non-linear evolution. As the gas flows 
have relatively low densities and moderate velocities, they 
can easily reach supersonic regimes. Therefore, they can oc- 
cupy an important volume of the simulated region. However, 
as it is shown in the right panel, these are weak shocks char- 
acterised by low volume- weighted Mach numbers. 

As the non-linear evolution of cosmic structures con- 
tinues, the number of collapsing regions increases, as well 
as the densities and velocities of the gas flows falling into 
such structures. The combination of these factors produces 
a maximum of the shocked volume and the shock strength 
around z ~ 5. 

After this time two different trends appear depending 
on whether we look at the shocked volume or at the average 
shock strength. The first one continuously declines mainly 
due to the fact that structures start to collapse and to form 
bound haloes where the gas component is locked up and, 
therefore, the sizes of regions where shocks can appear get 
smaller as the time evolution goes on. On the other hand, 
the average shock strength shows the same behaviour imme- 
diately after its maximum but it changes its derivate around 
z ~ 2. The reason of this behaviour is strongly correlated 
with the increase of the merging rate among the already 
formed haloes of different masses and sizes. The final fraction 
of the simulated shocked volume reaches a value of ~ 20% 
at 2 = 0. 
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Figure 1. Large-scale distribution of several quantities at 2 = 0. Each two-dimensional panel is 64 Mpc side length. All quantities are 
integrated along a slice of thickness L = 10 Mpc according with gi = / qdx where tj is a generic quantity to be projected and is 
the final plotted variable. The panels show: the integrated gas overdensity, p/pg, being p^ the background density (upper left), the 
integrated gas temperature in K (upper right), the integrated gas entropy in keV cm? (lower left) and the integrated Mach number (lower 
right). Different palettes have been used in order to highlight the particular features of each panel. 



This result on the temporal evolution of the fraction 
of the computational volume occupied by shocked cells 
se ems to be in broad agreement with the results shown 
in lVazza. Brunetti fc Ghelled (|2009l ) who obtained that the 
fraction of the simulated volume being shocked goes from 
roughly 30% at high redshift [z ^ 6) down to ~ 15% at 
z = Q. 

A crucial tool in the study of cosmic structures is the 
mass distribution function of the different objects in the 
Universe. Following a similar approach to the mass func- 
tion, we study the shocks distribution function (from now 
on SDF), defined as the variation of the number of shocked 
cells (volume-weighted) as a function of the Mach number. 
In Fig. 13] the differential SDF at 2 = is plotted - in log- 
arithmic scale - to show the distribution of shocks in the 
simulated volume. 

The most remarkable feature of the SDF is the presence 
of two clearly differentiated regimes with different slopes. 
The first regime corresponds to low Mach numbers which are 
much more abundant. The second regime is steeper than the 
previous one and stands for strong shocks with higher Mach 
numbers which are the most rare in the studied volume. 



The transition between the two regimes takes place around 
M r-u 20. This value for A4, as it will be discussed in the 
next sections, turns out to be a crucial number separating 
the internal shocks - located within the virial radius of the 
structures - and the external shocks at outer regions. The 
two different regions of the SDF can be fitted by power laws 
of the form dN{M)/dM oc Thus, for low Mach num- 

bers (up to ~ 20) we obtain a slope of a ~ —1.7, whereas 
for stronger shocks a steeper relation, a ~ —4.1, is found. 
Previous works already proved that the bulk of shocks in the 
universe is dominated by relatively weak shocks {M ^ 2), 
but there is also a population of stronger shocks located 
in the external regions of haloes where structures are not 



completely virialised (e.g..lRvu et al.l 20031 : iPfrommer et al.l 



completely vmaliseg (e.g..lKvu et al.l a 
12009 : IVazza. Brunetti fc GhellejbOOa ) 



We have also followed the time evolution of the SDF 
at several redshifts, finding very similar trends at all the 
epochs. The most relevant difference between the distribu- 
tions obtained at different times is found at the high-Mach 
end of the relation, which slowly moves down to lower values 
when the redshift decreases. Only when quite high redshifts 
(z ^ 6) are reached, an important change in the shape of 
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Figure 3. Shocks distribution function (SDF) at 2 = as a 
function of their Mach number in logarithmic bins. 

the distribution is appreciated although, probably, the in- 
formation provided by the SDF at such high redshift is not 
relevant due to the early stage of evolution of the cosmic 
structures. 

3.2 Shock waves and cosmic structures 

Cosmological shocks derived from the formation and evolu- 
tion of cosmic structures have to present some correlations 
with the main features of the distribution of such structures. 
In order to deepen in this connection, in the present Section 
we discuss the correlation between the distribution of cos- 



mological shock waves and the halo population (galaxies and 
galaxy clusters). 

The first analysis one can think about is to compare 
the shock pattern with the density distribution of the differ- 
ent components forming the haloes in the simulation. Thus, 
Fig. [4] shows a 2-D projection along the z axis of the Mach 
number distribution (upper left panel) to compare with the 
gas (p/Pb), dark matter {poM / p^) ^-i^d stellar [pt/p^) over- 
densities (upper right, lower left, and lower right panels, re- 
spectively) at 2; = 0. Each panel represents a thin slice of 0.2 
Mpc thickness and 64 Mpc side length centred at the posi- 
tion of the most massive halo found in the computational 
box (the same than in Fig.[T]). All the plotted quantities are 
in logarithmic scale. In all these panels, the contours of the 
strong shock waves - with high Mach numbers {M > 20) - 
are overplotted. 

From Fig. |4] we can infer, as it was expected, that 
the shock pattern perfectly traces the cosmic web (e.g., 
iMiniati et aDl2000l l. Nevertheless, this shock pattern has 
two clearly separated regimes. The first of these regimes 
gathers the high M shocks - overplotted as contour lines 
in Fig. [4]- that wrap filaments, sheets, and haloes. As these 
shocks are located out of the virial radius of the structures, 
they are classified as external shocks. These shocks have 
quasi-spherical geometries and can be located at distances 
of several virial radius f rom the centre of the structure w here 
they were created (e.g.. lVazza. Brunetti fc Ghellerll2009l '). In 
general, such external shocks are the evolved state of the ac- 
cretion shocks formed during the collapse of haloes. During 
their movement outwards from the halo centre, these accre- 
tion shocks interact among them creating, therefore, a more 
complex pattern. 

Although for the sake of clarity the low-Mach number 
shocks are not plotted in Fig. |4l moving inwards the virial 
radius of haloes, more irregular and weaker shocks {M ^ 5) 
are formed, in line with results from previous studies (e.g., 
IVazza. Brunetti" fc Gheller|[2009l ). This second regime corre- 



8 S. Planelles & V. Quilis 




Figure 4. Distributions of Mach number compared with dark matter {pDM I P b ) > {p/ Pb ) ^^'^ stellar {pt/ ) overdensities at 2 = 0. 
Each panel is a thin slice of 0.2 Mpc thickness and 64 Mpc side length. They show the Mach number distribution (upper left) and the 
gas, dark matter and stellar overdensities (upper right, lower left, and lower right panels, respectively). In all the panels, the contours of 
the shock waves with high Mach numbers are overplotted. 



sponds to the internal shocks which are mainly associated to 
random flow motions and merger events within the haloes. 

In order to reinforce the previous idea, in Fig. [5] we 
compare the distribution of shock waves with high Mach 
numbers (given by the contour lines) with the spatial distri- 
bution of the dark matter haloes found in the simulation at 
z = Q. This figure shows the same slice than in Fig. |4l The 
represented dark ma tter haloes have been iden tified with the 
ASOHF halo finder (|Planelles fc Quilidl201ol ). Only objects 
with virial massed larger than 10^^ Mq are considered. The 
total sample consists of ~ 260 haloes at 2 ~ within a 
range of masses between 1.0 x 10^^ Mq and 8.0 x lO" Mq. 
The dark matter haloes are plotted as circles whose sizes 
represent their virial radius. Those haloes whose position 
perfectly fits in the analysed slice are drawn in red. It is 
clearly visible how the complex pattern of strong shocks - 



The virial mass of a halo, M^ir, is defined as the mass en- 
closed in a spherical region of radius with an average density 
Ac = ISvr^ -I- 82x — 39x^ times the critical density of the Universe, 
being x = 0(2) - 1 and ^(z) = [^^(1 + zf]/lnm(l + z}'-^ + ^a] 
llBrvan fc Norman|[l998h . 



produced by the hierarchical evolution of the proto-haloes 
- surrounds the galaxy cluster (large red circle) from which 
these shocks stem from. 

3.3 The Mach number scaling relation 

The scaling relations are powerful tools in the understand- 
ing of the processes sculpting the formation and evolution 
of galaxy clusters and groups. The common way to use the 
scaling relations consists in plotting, for a given character- 
istic radius, the cluster mass against the temperature, the 
entropy, or the X-ray luminosity. 

Starting from this idea, we generalize its fundamental 
concept and we introduce what we call the Mach number 
scaling relation. In order to build this relation we compute 
an average volume-weighted Mach number for each halo in 
our sample with Af„i,. > IO'^^Mq. This volume- weighted 
Mach number is computed within the virial radius of each 
halo and, therefore, it represents a virial Mach number. After 
that, we bin individual haloes into a two-dimensional grid of 
Mach number versus virial mass in order to have information 
on the number of haloes within a certain range of Mach 
numbers and viral masses. The obtained scaling relation is 
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Figure 6. Halo distribution in the Mach number - virial mass plane. Results at z ~ 2, 1, 0.4, and 0, are shown. The shaded regions have 
been computed by binning individual haloes into a two-dimensional grid in the plane M — Myi^. Six contour lines equally spaced are 
plotted to highlight the two-dimensional distributions. 



displayed in Fig. [6] in which we show different contour lines 
in order to highlight the shape of the 2-D distribution. Four 
different epochs corresponding to z ~ 2,1,0.4, and 0, are 
shown. The sample of haloes with Mu,> > 10^^ Mq at these 
epochs ranges from ~ 160 haloes at 2; ~ 2 to ~ 260 at 2 = 0. 

The analysis of Fig. [6] revels some striking features. Let 
us focus first on the analysis of the distribution at z = 
(lower right panel), where there are two well-differentiated 
trends. On the one hand, there is an almost constant region 
for low Mach numbers (Ai < 5) which seems to be no- 
dependant on the mass of the haloes. On the other hand, a 
steeper trend along all the range of Mach numbers seems to 
be correlated with halo masses. 

Our explanation for this distribution is that the first 
trend is occupied by haloes which started their evolution 



in a relatively smooth way, with no mergers or a few mi- 
nor mergers. The shocks within the virial radius producing 
the average Mach numbers associated with these haloes are 
related with smooth accretion flows of gas falling into the 
objects during this quiescence evolution. Therefore, haloes 
that have been quietly set up since early phases of their evo- 
lution have low Ml an tend to be located at the bottom of 
the plane M — Mvir independently of their virial masses. 

Complementary, haloes that initially were involved in 
merger events and in active phases of their evolution suffer 
violent events that produce stronger shocks compared with 
those in the previous region. In this branch of the M — Mvir 
plot, there is a strong dependence on the virial masses of 
haloes. 

Therefore, the evolutionary history of each halo explains 
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Figure 5. Distribution of dark matter haloes together with shock 
waves at z = 0. The sohd Hnes represent six contour levels of the 
higher-Mach number shock waves. This contour has been com- 
puted for a slice of the simulated box of 0.2 Mpc depth and 64 
Mpc side length centred at the location of the most massive halo. 
Circles stand for all the found haloes in the computational box 
with masses larger than ^ 10^^ A^q. We point out in red those 
haloes whose positions fit in the considered slice. The size of the 
circles represent the virial radius of the haloes. 



the bimodal segregation shown in Fig. [6] which is also ob- 
served through the temporal evolution. At early epochs, cor- 
responding to the formation time of large galaxies and small 
groups of galaxies {z ^ 2), an L-like pattern in the plane 
M — Mvir appears. These times correspond to highly non- 
linear epochs of the evolution of the haloes, with high rates 
of merger events and interaction among the structures at dif- 
ferent spatial and mass scales. When advancing in time, the 
initial L-pattern trend progressively bends over to higher 
masses - corresponding to big galactic haloes and galaxy 
clusters - but lower Mach numbers, reaching the bimodal 
distribution, previously discussed, at 2 ~ 0. 

In order to understand this behaviour, we have followed 
the global evolution of individual haloes looking at their 
overall drifts along the plane M — M^jir- We have studied 
this evolution for the 22 more massive haloes in the simula- 
tion which, indeed, are the best numerically resolved. The 
haloes in this subsample evolve according to two different 
behaviours. Roughly the 50% of this subsample of haloes 
begin, at high redshift, with a relatively high Mach num- 
ber and evolve to progressively lower Mach numbers while 
increasing their mass. The remaining percentage of haloes 
depart from low-Mach number states {M < 5) and tend to 
move, during their evolution, almost parallel to the x-axis 
while augmenting their masses. Since our sample of haloes 
is far from being statistically complete we can not make a 
robust conclusion. However, our hypothesis to explain this 
behaviour is that the evolution of haloes along the plane 
M — Mvir is intimately related with their dynamical his- 
tory. Like a gross trend, those haloes suffering important 



merger events (major mergers) early in their evolution only 
can evolve towards lower Mach numbers while reaching an 
equilibrium state producing, therefore, the decline of the ini- 
tial L-like pattern into the flatter final distribution. On the 
other hand, haloes with a relatively quiet beginning show 
very low initial Mach numbers and evolve without signifi- 
cant changes in their virial Mach number and, consequently, 
they move almost parallel to the M„ir-axis while increasing 
their mass. 

Let us point out that this general behaviour is not 
in contradiction with previous studies. Once haloes, due 
to their early evolution, are initially placed in the M — 
Mvir plane, shocks {A4 ^ 3) inside their virial radius 
can be episodically found associ ated to merger events (e.g., 
IVazza. Brunetti fc GhelleJ 120091 '). However, the strength of 
these internal shocks is not enough to break the shape of the 
obtained L-like pattern. 

Figure[B]also encloses relevant information on the evolu- 
tion of global integrated quantities. In this regard, the time 
evolution of the mean Mach number of the overall sample 
of haloes clearly shows a progressively decline from higher 
values {Ai — 10) at high redshifts, down to lower values 
{M ~ 5) at 2 = 0. This decrease of the mean Mach number 
is fully consistent with a picture of haloes evolving towards 
an equilibrium state. 

This behaviour can be correlated with the one shown in 
Fig. [3 In this figure, we show the evolution with redshift of 
the mean mass-weighted virial Mach number for the entire 
population of haloes at a given epoch. In this analysis we 
are dealing with internal shocks and, therefore, the Mach 
numbers are relatively low, going from — 25 at z ~ 6 
down to ~ 5 at 2 = 0. As when discussing Fig. [2l it is 
important to mention that resolution effects could affect the 
results, especially at high redshift. 

In any case, and independently of the considered epoch, 
the mean Mach number of the bulk of haloes, as derived from 
Figs.[6]and[7l is always lower than ~ 20 — 25. This result per- 
fectly correlates with the turn observed in the distribution 
function of shocks discussed in Sec. 3.1 (see Fig. (3]), repre- 
senting the transition between external and internal shocks. 

3.4 Halo properties and shock waves 

The results presented in Sec. 3.2 reveals the already ex- 
pected conclusion that there should exist some sort of cor- 
relation between the external shocks surrounding clusters 
and galaxies and the dark matter haloes where such shock 
waves were probably originated. In this Section, we are going 
to study this correlation in detail and its possible implica- 
tions to characterise the haloes where the shock waves are 
originated. 

The collapse of an isolated and idealised halo would 
produce an accretion shock that, after the formation of the 
halo, would move outwards from the centre of the structure 
defining an spherical shell. In a more realistic scenario, and 
assuming that the formation of every halo produces an ac- 
cretion shock, the shocks surrounding the haloes would not 
be perfectly spherical and in so me cases the formed patterns 
can be extremely complex (e.g.. iMiniati et al]|2000l ). Never- 
theless, there are situations where these shocks can be well 
defined even in a realistic case: i) for those haloes which have 
had a relatively quiet history, and ii) for the most massive 
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Figure 7. Mean mass- weighted virial Macii number for the entire 
population of haloes as a function of redshift. 

haloes which usually have stronger accretion shocks at their 
outskirts than the smaller objects. 

In Fig. [8] we display several images showing the 
distribution of shock waves at several redshifts (2 = 
1.4,1,0.5,0.2,0.13,0) around the most massive halo in the 
simulation. All the panels represent slices of 0.2 Mpc thick- 
ness and ~ 21 Mpc side length. In each of the panels in 
Fig. [SI in addition to the shock pattern, we overplot the 
position and size of the dark matter halo as a white circle 
whose radius represents the halo virial radius in comoving 
coordinates. At z = 0, the quasi-spherical accretion shock 
wave associated to the halo is perfectly visible and easily 
identifiable. By simply tracking backwards in time in the 
different panels the position of the shock, we can arrive at 
the time when this accretion shock was formed. If we assume 
that the formation of the accretion shock can be defined as 
the formation time of the halo, this procedure produces a 
new and natural way to estimate the formation time of a 
galaxy cluster, in contr ast to the common d efinition of the 
cluster formation time (|Lacev fc Cold 1 19931 ). i.e., when the 
halo reaches half of its mass at z = 0. 

Besides the qualitative procedure just described, we 
have tried to track the movement of the accretion shock 
around the haloes in greater detail. In order to do so, we have 
visually identified and measured the most relevant quanti- 
ties of the shocks. Unfortunately, the statistical limitation of 
our sample of massive haloes, plus the difficulty of the semi- 
manual methodology to follow the evolution of the accretion 
shocks, has made extremely difficult to perform this analysis 
and, as a consequence, we have been able to do it only for 
the most massive halo in the simulation. In Table [1] we sum- 
marise the main data corresponding to the studied halo and 
its associated external shock wave according to Fig. [8] The 
measured quantities are: the redshift (z), the cluster virial 
meiss (M„ir), the cluster virial radius {Rvir), the comoving 
radius of the shell defined by the shock {Rs), and the shock 
speed {vs). The radial distance from the halo centre to the 



Table 1. Main data describing the most massive halo and its 
associated external shock wave. The different columns stand for 
the redshift, z, the halo virial mass, M^ir, in units of 10^* Af0, 
the halo virial radius, R^ir, in comoving Mpc, the radius of the 
shock wave, Rs, in comoving Mpc, and the shock speed, Vs, in 
km/s, respectively. 
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shock, Rs , is estimated using a visual mechanism which can 
only produced approximated results. To do so, we make thin 
slices of the Mach number distribution like the ones plotted 
in Fig. [S] In these slices, we identify the halo and associate 
by eye the thin quasi-spherical shell having high and ho- 
mogenous Mach numbers to the shock wave that wraps the 
halo. Then, considering several directions, we measure dif- 
ferent values for the radial distance from the centre of the 
halo to the shock front. The final value for the shock ra- 
dial distance, Rs , is obtained as the arithmetic mean. Using 
the same procedure, the modules of the gas velocity at the 
different points of the shock shell used to estimate Rs, are 
computed and averaged to produce the final shock speed. 

Vs. 

The propagation of the shock wave surrounding the con- 
sidered halo can be well summarised in Fig.[51 where the dis- 
tance from the shock wave to the halo centre {Rs) is plotted 
against the time. As expected, the behaviour is consistent 
with the self-similar nature of an spherical isolated shock 
wave, being the slope of the line plotted in Fig. [9] the prop- 
agation speed of the shock. 

The last column in Table [T] shows the shock speed {vs) 
at different times. These values are consistent with a con- 
stant propagation speed. Nevertheless, there are deviations 
from a perfectly constant slope which are explained by the 
errors associated to the graphical method used to measure 
the positions and velocities of the shock wave, and by the 
non idealised environment surrounding the haloes where the 
external medium is far from being homogeneous. 

Assuming all these uncertainties, we can use Fig. [9] to 
date the formation time of the studied halo around z ~ 1.4 
and, therefore, it is nowadays ~ 8.9 Gyrs old. This redshift is 
the highest one at which the shock appears perfectly visible 
and well defined. At early stages, the situation is too messy 
and the crude method used to track the shock produces too 
much error. In this sense, this time should be taken as a 
rough estimate of the collapse of the halo. In any case, it is 
interesting to compare this formation time of the halo with 
its half mass formation time which for this particular cluster 
is 2 ~ 0.2. 
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Figure 9. Temporal evolution of the shock wave distance to the 
halo centre. 

Under the hypothesis that the shocks associated with 
haloes could be observed and measured, we suggest the idea 
to use them as probes to infer the main characteristics of 
the haloes where they were created. Although the detection 



and measurement of large-scale shock waves is an observa- 
tional challenge, there is hope that new com ing instruments 
- i.e. , SKA Q (Square Kilom etre Array; e.g.. iRawlings et al.l 
12004 iBattagUa et al.ll2009^ - could shed some light on this 
scenario. 

As it has been already mentioned, the properties of the 
shock waves have been derived by a graphical method. Be- 
sides, our sample is very reduced. Therefore, due to the un- 
certainties of the method and the statistical limitations, the 
results of this Section should be interpreted with caution, 
and only as a way to illustrate the potentiality of the study 
and observation of the shock waves associated to cosmic 
structures. 



4 SUMMARY AND CONCLUSIONS 

In this paper, we focus on the properties of the shock waves 
developed during the hierarchical evolution of the cosmic 
structures by means of a cosmological AMR simulation. 

To analyse the shock waves, we develop a numerical al- 
gorithm able of identifying shocks in 3-D AMR simulations. 
After labeling all the shocked cells (compression regions with 
V ■ V < 0) within the computational box, the Mach num- 
bers of the shocks are computed by means of the Rankine- 
Hugoniot temperature-jump condition. The Eulerian nature 
of the adopted numerical scheme allows to accurately detect 

http://www.skatelescope.org 
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shocks. Besides, the AMR structure of the employed numeri- 
cal method produces sets of grids mapping different levels of 
resolution that naturally permits to find shocks correspond- 
ing to different spatial scales. 

We analyse the morphology of the shock patterns 
detected in the simulated volume which turns out to 
be rather comple x. In qual i tative agreement with previ- 
ous works (e.g.. IRvu ct al.' '2Q03|; IPfrommer et al.l bood : 
IVazza. Brunetti &l Ghellcr. ,20091 ). the shocks morphologies 
follow the shape of the cosmic web, being filaments, sheets, 
and haloes surrounded by strong external shocks, while more 
irregular and weaker internal shocks are found within the 
virial radius of haloes. From the morphological point of 
view, it is remarkable that the external shocks associated 
to the most massive haloes - galaxy clusters - show a quasi- 
spherical shape around them. 

The statistical analysis of the shocks found in the simu- 
lation has produced several important conclusions. We have 
estimated that the filling factor of shocks - the fraction of the 
total volume occupied by shocked cells - at 2 = is roughly 
the 20% of the simulated volume with a mean Mach number 
of ~ 4. 

Following with this analysis, we have computed the dif- 
ferential shock distribution function within the simulated 
volume as a function of the shock Mach numbers. This func- 
tion can be fitted by two different power laws with the 
form dN{M)/dM oc TVf™, having the low Mach numbers 
(up to ~ 20) a slope of a ~ —1.7, whereas a steeper re- 
lation (a ~ —4.1) stands for stronger shocks. This turn 
has a crucial meaning as it shows the transition between 
the scales associated to internal and external shocks. The 
distribution function also reveals that most part of cosmo- 
logical shocks are essentially weak shocks (M ^ 2), as al 



ready demons trated b y other authors (e.g., Ryu ct al."2003l: 
IPfrommer eTa l. 2006: ,Vazza. Brunetti fc Ghellcr 2009). 

In agreement with the information derived from the 
shock distribution function, we find that the average Mach 
number within the virial radius of haloes at 2 = is ~ 5. 
If we look at higher redshifts, this average Mach number is 
always below ~ 20 — 25, which perfectly correlates with the 
turn obtained in the shock distribution function. 

When haloes, according to their properties, are located 
in the plane formed by the two axes, virial mass and aver- 
age Mach number inside the virial radius, we observe two 
well-separated trends at 2 = 0. One stands for a group of 
haloes with almost constant low Mach numbers (up to 5) 
which seems to be non dependent on the halo mass. This 
population of haloes, characterised by low virial Mach num- 
bers, could stem from the primordial formation of struc- 
tures which have been initially set up in a quasi-relaxed and 
smooth process. 

The second group of haloes spreads all over the range 
of Mach numbers but showing a clear correlation with the 
halo masses. They represent the haloes where strong shocks 
took place during their early evolution as a consequence of 
early mergers and violent accretion processes. We conclude 
that the evolution of haloes along the plane A4 — M„ir is 
intimately related with their merging history and with their 
initial set up. In that sense, possible estimates of Ai for a 
given halo could lead to know whether it has suffered merger 
events in the past, and it may correlate with other relevant 
features like the existence or not of cool cores. 



Finally, we speculate about the possibility that new ra- 
dio observations may observe the shocks associated with the 
formation and evolution of the cosmic structures. Assuming 
that this could be feasible in the near future, we propose to 
turn the argument around, and to use the main features of 
shocks - position and strength - to infer import features of 
the halo harbouring the shock like, for instance, the forma- 
tion time of the halo, or its total mass. This last application 
could be a powerful new tool for cluster mass estimation if 
we are able t o detect the accretion s hocks with radio obser- 
vations (e.g.. iGiacintucci et akllioOSl ). 

Our main conclusion is that shock waves play a crucial 
role in the formation and evolution of galaxies and galaxy 
clusters. Despite this general and obvious statement, direct 
evidence of shocks, both from the cosmic web formation pro- 
cesses (large scales) and those due to cluster merging events 
(small scales), has been found only in a relatively small num- 
ber of clusters thanks to observations of radio relics and 
temperature maps in X-rays. Therefore, we believe that it 
is necessary to pursue the study of the role of shocks in 
cosmological context as they are fundamental players in the 
paradigm of structure formation in the Universe. Their com- 
plete theoretical description together with their detection 
and observation are still a challenge. 
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